! program flap4rcr
!     use mod_param
!     implicit none
    
!     real :: th0, th1c, th1s
!     real :: beta0, beta1c, beta1s 

!     th0 = 0.14
!     th1c = -0.06
!     th1s = -0.04

!     call sol_beta(th0, th1c, th1s, beta0, beta1c, beta1s)

!     print *, "The 1st flap response is:", beta0/pi*180.0, "+", &
!              beta1c/pi*180.0, "cos(Phi)", &
!              beta1s/pi*180.0, "sin(Phi)"
    
!     read(*,*) 

! end program flap4rcr

subroutine sol_k10 (K10)
    use mod_param
    implicit none 

    real :: K10 

    ! Omega = 35
    ! omegan = 1.4
    ! e = 0.25 
    ! m_ = 8.114
    ! Ms = 0.5*m_*R**2.0
    ! Ib = 1.0/3*m_*R**3.0

    K10 = (omegan**2-1-(e*Ms)/Ib)*Ib*Omega**2

    return 
end subroutine sol_k10

subroutine sol_v1d(thtot, lam1, v1d)
    use mod_param
    implicit none 

    real :: thtot ! input 
    real :: lam1, v1d ! output 
    real :: ctt, lamf ! tmp
    ! real :: v1d0 ! tmp

    ! vair = vin
    ! T = 5500.0

    ! Omega = 35.0
    ! R = 5.49
    ! chord = 0.365
    ! sigma = 0.0635
    ! Arot = pi*R**2
    ! cl_a = 6.0
    ! as = -3.0/180*pi 

    ! vtip = Omega*R 
    ! CT = T/2.0/(0.5*rho*vtip**2*Arot)
    ! mu = vair*cos(as)/vtip
    ! lam0 = vair*sin(as)/vtip

    ! v1d0 = CT*vtip/(lam0**2+mu**2)
    ! lam1 = lam0-v1d0/vtip
    ! do while(.true.)
    !     v1d =  CT*vtip/sqrt(lam1**2+mu**2)
    !     if(abs(v1d-v1d0)<=1e-3) then
    !         exit 
    !     else
    !         v1d0 = v1d
    !         lam1 = lam0-v1d/vtip
    !     end if
    ! end do 

    if(abs(vair-0.5)<=1) then
        call sol_ct4vf(thtot, ctt, lamf)
        print *, "v1d VF: ", lamf
    else 
        call sol_ct4ff(thtot, ctt, lamf)
        print *, "v1d FF: ", lamf
    end if 

    lam1 = lamf 
    v1d = (-lam0+lamf)*vtip

    return 
end subroutine sol_v1d

subroutine sol_beta(th, beta)
    use mod_param
    implicit none

    real :: th(1:3) ! 输入参数
    real :: beta(1:3) ! 返回参数
    real :: K10, v1d, thtot, &
            ctt, lamf, lam1 ! 内部使用
    
    ! T = 5500

    ! Omega = 35.0
    ! R = 5.49
    ! chord = 0.365
    ! sigma = 0.0635
    ! Arot = pi*R**2
    ! cl_a = 6.0
    ! as = -3.0/180*pi 

    ! vtip = Omega*R 
    ! CT = T/(0.5*rho*vtip**2*Arot)
    ! mu = vair*cos(as)/vtip
    ! ! lam0 = vair*sin(as)/vtip

    ! eta0 = 7.5/180*pi 
    ! thtw = -10.0/180*pi 
    thtot = th(1)

    ! m_ = 8.114
    ! Ms = 0.5*m_*R**2.0
    ! Ib = 1.0/3*m_*R**3.0
    ! lockb = 5.41
    ! omegan = 1.4

    call sol_K10(K10)
    call sol_v1d(thtot, lam1, v1d)
    print *,"v1d is calculated!"
    if(abs(vair-0.5)<=1) then
        call sol_ct4vf(thtot, ctt, lamf)
        print *, "VF is calculated!"
    else 
        call sol_ct4ff(thtot, ctt, lamf)
        print *, "FF is calculated!"
    end if 

    E = (144/lockb**2*(e/R)**2)/((1-e/R)**6*(1+0.5*mu**2))

    beta(1) = (((3.0/8)*lockb*(ctt/sigma*cl_a)*(1-e/R)**2)/(1+0.5*e/R+K10/Omega**2)) &
            -((1.5*g*R/vtip**2)/(1+0.5*e/R+K10/Omega**2))
    beta(2) = (2.0*mu*(4.0/3*thtot+thtw+mu*as-v1d/vtip)-(1.0+1.5*mu**2)*th(3)) &
            /((1.0-0.5*mu**2)+E) &
            + (((12/lockb*e/R)*(th(2)*(1+0.5*mu**2)+8.0/9*(ctt/(cl_a*sigma))*(lockb*mu/(1+1.5*e/R)) &
            + 4.0/3*v1d/vtip))/((1-e/R)**3*(1+0.5*mu**2))) &
            /(((1-0.5*mu**2)*E))
    beta(3) = (1.0/(1+0.5*mu**2)*((4.0/3*mu*beta(1)/(1+1.5*e/R))+v1d/vtip)+th(2)) &
            -(12.0/lockb*e/R*(2*mu*(4.0/3*thtot+thtw+mu*as-v1d/vtip)-(1+1.5*mu**2)*th(3))) &
            /((1-e/R)**3*(1-0.25*mu**4))

    return
end subroutine sol_beta

subroutine sol_bint(th0, th, beta)
    use mod_param
    implicit none

    real :: th0(1:3), th(1:3) ! 输入参数
    real :: beta(1:3) ! 返回参数
    real :: K10, v1d, thtot0, thtot, &
            ctint, lamint, &
            lam10, v1d0! 内部使用
    
    ! vair = vin

    thtot0 = th0(1)
    thtot = th(1)

    call sol_K10(K10)
    call sol_v1d(thtot0, lam10, v1d0)
    call sol_ct4int(thtot0, thtot, ctint, lamint)
    print *, "CT for Intrupt is calculated!"

    E = (144/lockb**2*(e/R)**2)/((1-e/R)**6*(1+0.5*mu**2))

    beta(1) = (((3.0/8)*lockb*(ctint/sigma*cl_a)*(1-e/R)**2)/(1+0.5*e/R+K10/Omega**2)) &
            -((1.5*g*R/vtip**2)/(1+0.5*e/R+K10/Omega**2))
    beta(2) = (2.0*mu*(4.0/3*thtot+thtw+mu*as-v1d/vtip)-(1.0+1.5*mu**2)*th(3)) &
            /((1.0-0.5*mu**2)+E) &
            + (((12/lockb*e/R)*(th(2)*(1+0.5*mu**2)+8.0/9*(ctint/(cl_a*sigma))*(lockb*mu/(1+1.5*e/R)) &
            + 4.0/3*v1d/vtip))/((1-e/R)**3*(1+0.5*mu**2))) &
            /(((1-0.5*mu**2)*E))
    beta(3) = (1.0/(1+0.5*mu**2)*((4.0/3*mu*beta(1)/(1+1.5*e/R))+v1d/vtip)+th(2)) &
            -(12.0/lockb*e/R*(2*mu*(4.0/3*thtot+thtw+mu*as-v1d/vtip)-(1+1.5*mu**2)*th(3))) &
            /((1-e/R)**3*(1-0.25*mu**4))

    return
end subroutine sol_bint 